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Abstract 

We present a model for mechanically-induced pattern formation in 
growing biological tissues and discuss its application to the development 
of leaf venation networks. Drawing an analogy with phase transitions in 
solids, we use a phase field method to describe the transition between 
two states of the tissue, e.g. the differentiation of leaf veins, and consider 
a layered system where mechanical stresses are generated by differential 
growth. We present analytical and numerical results for one-dimensional 
systems, showing that a combination of growth and irreversibility gives 
rise to hierarchical patterns. Two-dimensional simulations suggest that 
such a mechanism could account for the hierarchical, reticulate structure 
of leaf venation networks, yet point to the need for a more detailed treat- 
ment of the coupling between growth and mechanical stresses. 

1 Introduction 

A broad range of patterns in nature are induced by mechanical forces. While 
the most well understood examples are found in physical phenomena, e.g. in 
solids, fracture [T], wrinkling [2], delamination [3 , or localized deformation [5], 
it is increasingly becoming appreciated that mechanical forces are also essential 
in the morphogenesis of living organisms. In the simplest case, biological tissues 
can undergo mechanical instabilities as a result of inhomogeneous growth [5l [6] . 
More broadly, the importance of the mechanical and adhesion properties of cells 
in morphogenesis [7] and the ability of mechanical forces to act as a signal in 
growing tissues and organisms [5J EE HO] suggest that the interplay between 
mechanical forces and biochemical factors has a central role in development. 

An illustration of these ideas can be found in recent investigations of the 
structure of leaf venation networks [llj . The standard model for vein forma- 
tion involves the transport of the plant hormone auxin |12) : auxin circulates 
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through the leaf, with the leaf base acting as a sink; as the leaf develops, auxin 
flow becomes progressively concentrated along preferential paths, which mark 
the location of future veins. This model has received ample experimental sup- 
port (see e.g. [13]), yet it fails to explain certain aspects of leaf venation, such 
as the presence of closed loops. Simulations of auxin transport generally give 
rise to tree- like structures |14] , a typical property of diffusion- limited growth 
processes [15]. In contrast, it has been noted that the reticulate structure of 
venation networks is very similar to that of crack patterns, suggesting the pos- 
sibility that the differentiation of leaf veins could be governed by mechanical 
stresses [TT]. The inner tissues of the leaf, where the veins form, are subjected 
to compressive stresses by the epidermis, and one of the earliest signs of vein 
differentiation is the elongation of vascular cells [16], which could release the 
compressive stresses in their vicinity, providing a feedback between differentia- 
tion and the mechanical state of the tissue. 

In this article, we develop a theoretical approach to study the formation 
of patterns induced by mechanical stresses in biological tissues. Of particular 
interest is the fact that these tissues are growing as the patterns form. We show 
how growth can be the driving force that causes the formation of new structures, 
and how the history of their formation can be reflected in their organization. 
Let us emphasize that this study mainly aims at demonstrating possible effects 
of mechanical stresses, and that the approach developed here would have to be 
integrated with biochemical mechanisms in a more complete model. 

Stated in general terms, the patterning mechanism we consider is the follow- 
ing: compressive stresses imparted on a tissue cause it to switch to a "collapsed" 
state occupying a smaller volume, whereby the compression is released. In the 
case of leaf venation, this transition describes the differentiation of veins, and 
the collapsed state refers to the elongated vascular cells. An analogy can be 
made with the behavior of solid foams, which exhibit localized deformation 
upon compression [4] . Our model is not explicitly constructed according to this 
analogy. Instead, we consider the initial and collapsed states of the tissue as two 
phases of a continuous medium, which we describe using a phase field method. 
However, the response of this two-phase medium is analogous to that of a non- 
linear elastic material, and can be analyzed with reference to the behavior of 
such materials 0. 

In what follows, we first present our model in a one-dimensional setting. 
We define a minimal model, amenable to analytical treatment, that implements 
the patterning mechanism described above. We consider a layered system, in 
which the tissue is coupled to a rigid substrate, which could represent a stiffer 
tissue such as the epidermis of plant leaves. Stresses are induced by a mismatch 
between the two layers, which can occur as a result of differential growth. If 
the transition between the two states of the tissue is reversible, regular patterns 
are obtained, in which domains of collapsed and uncollapsed tissue alternate. 
If, instead, the transition to the collapsed state is made irreversible to model 

1 Conversely, localized deformation in non-linear elastic materials can be interpreted as a 
first-order phase transition. 
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tissue differentiation, hierarchical patterns are obtained. 

When turning to two-dimensional systems, the tensorial nature of mechani- 
cal fields comes into play, and a much greater variety of behaviors are possible. 
Rather than to attempt a systematic exploration, we present a few illustrative 
examples of the patterns that can be obtained. In particular, we show that 
anisotropic collapse, which could describe the elongation of vascular cells in 
plant leaves, produces reticulate patterns comprised of interconnected stripes 
of collapsed tissue. Our simulation results suggest that hierarchical, reticulate 
patterns might be obtained when irreversibility is introduced. However, in the 
simple model presented here, these patterns are disrupted by strong residual 
stresses that develop in the vicinity of growing collapsed regions. We suggest 
that a more satisfactory model would require a more detailed description of the 
interplay between mechanical stresses and growth. We also contrast our model 
with a similar approach presented in [17j . 

2 One-dimensional model 

In this section, we present our modeling approach in the context of a one- 
dimensional system. We consider a tissue that can exist in two states, an "ini- 
tial" state and a "collapsed" state that has a smaller volume. Transition to the 
collapsed state is triggered by compressive stresses, and releases the compres- 
sion. While, in general, the two states of the tissue could differ in their stiffness 
(the collapse could be due to the transition to a softer state), we assume that 
they only differ in their rest configurations (the collapsed state has a smaller 
volume at rest). In the case of leaf veins, the initial and collapsed states of the 
tissue represent the undifferentiated and vascular tissues, respectively, and the 
change in rest configuration corresponds to the elongation of vascular cells as 
they differentiate. The tissue is described as a continuous medium, and its two 
states as two "phases" , using a phase field model. 

The phase field approach, used here in the context of mechanically induced 
pattern formation [17], was first introduced to describe the growth of a solid 
in a liquid [H [H [HI [5T] (reviewed in [H] H3J) and has later proved to be 
a powerful tool to describe free boundary problems in fluid flows [24] [25] and 
solid-solid phase transitions in alloys [26], [27] . This is especially the case in 3D, 
where interface tracking methods are extremely difficult to implement [55] [55] . 
Our model is similar to those developed for fracture in [30] [311 13H 1331 GH] and for 
solid-solid phase transitions in [2l)l 12"?] , which comprise phases having different 
mechanical properties. 

A continuously varying auxiliary field (f> is introduced to indicate the local 
state of the tissue, with the convention that <j> — corresponds to the initial 
state, and <j) = 1 to collapsed regions. The behavior of the system, i.e. both 
its deformations and the transitions between the two states, is assumed to be 
governed by an energy functional E{(j), u}, where u is the displacement field. 
While this is natural for a physical system, it may seem rather arbitrary for a 
biological one. However, we will see that the resulting behavior is sufficiently 
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general. The evolution of the system is assumed to be quasi-static, i.e. at any 
given time it is in a state of equilibrium that minimizes E, and satisfies the 
equilibrium equations 
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For a one-dimensional system, the energy E is defined by 



E = 



^(a^) 2 + e(<M + /W>) 



dx. 



(3) 



The first term in the integral is a squared gradient term, which penalizes sharp 
variations in the phase field, e is the density of elastic energy, which depends 
on the phase and on the strain e = d x u. f(<f>) is a potential that, together 
with the coupling between <f> and e through e, governs the behavior of the phase 
field. Assuming that both phases (treated as linear elastic) have the same elastic 
modulus fi and differ only in their rest configuration, the energy defined by eq.[3] 
takes the form 
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(d x 4>) 2 + £ (e - e r (0)) 2 + f{fl>) 



dx, 



where the rest configuration is given by 



e r {4>) = h{4>)e r i. 



(4) 
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The function h (defined below) satisfies h(0) = and h(l) — 1, and e r \ is the 
rest configuration of the collapsed phase. The potential / is defined by 



f{(f>) = ahtf) 2 + /3h(<t>). 



(6) 



h(4>) 



3(f) 2 — 2c/) 3 ' is chosen such that h'(Q) = h'(l) = 0, so that uniform phases 
((f) — or <f> — 1) are stationary states. 

Qualitatively, it can seen be from eq. 2] that under sufficient compression, 
the system will tend to switch to the collapsed state to lower its energy. The 
response of the system can be described by examining its uniform states of 
equilibrium (where <f> and e are uniform). In the range of parameters considered 
herein, (f> has a unique equilibrium value for every value of the strain e. Defining 
the critical strains 
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The parameters are chosen such that the critical strains defined below satisfy t c l < fcO- 
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Figure 1: (a) Sketch of the stress-strain curve of the model. The dashed segment 
corresponds to unstable states, which phase-separate. The dotted horizontal line 
indicates the plateau stress observed when the two phases coexist, (b) Stress- 
strain response of a simulated inhomogeneous system. 

4> is equal to when e > e c o (note that since we are considering compressive 
strains, e and a are generally negative), varies gradually from to 1 as e goes 
from e c o to e c i, and is equal to 1 when e < e c i (for strong compressive strains). 
Since <f> is uniquely determined by the strain e, so is the stress a = fx(e — e r (<£)). 
Our two-phase model thus behaves as a non-linear elastic material with a certain 
stress-strain response. 

As shown on fig. QJa), this response is non-monotonic, and the phase transi- 
tions in our model are equivalent to mechanical instabilities in an elastic material 
with a non-monotonic stress-strain curve [35 . The states for which the slope 
is negative, which correspond to values of the phase intermediate between 
and 1, are unstable. If the system is subjected to increasing strain, it deforms 
uniformly until e = e c o, then phase-separates. When the two phases coexist, 
the stress has a fixed "plateau" value a p , as do the strains in the two phases 
(sec fig. H{a)). In a biological context, the plateau stress can be understood as 
a "homeostatic stress" restored by the formation of a collapsed region. If the 
strain is further increased, the proportion of the collapsed phase increases until 
the system is entirely collapsed, after which the strain is uniform again. Note 
that states intermediate between the plateau stress and the instability thresh- 
old are metastable, and in actual physical systems, the distinction between the 
two can be blurred by disorder 0J. This can be reproduced in our model by 
adding a random prefactor r(x) to the terms e(</>, e) + /(</>) in the expression of 
the energy (eq. [3J. This creates local weak spots, so that the collapsed phase 
nucleates before the critical stress for an ideal, homogeneous system (compare 



As illustrated by the simple, piecewise linear stress-strain curve shown in 
fig. [Ha), our model provides a minimal description of a mechanism for the re- 
lease of mechanical stresses by localized deformation. There are just enough 
parameters to allow independent adjustment of the properties that will deter- 
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Figure 2: Diagram of the system. A layer of tissue (top) is coupled to a rigid 
substrate (bottom). 



mine its pattern-forming behavior, i.e. the instability threshold, the amount 
of deformation associated with the collapse, and the equilibrium stress at the 
interface between the two phases. 

In what follows, we consider the system shown in fig. [2J in which a layer of 
tissue described by the above model is coupled to a rigid substrate. If there 
is any mismatch between the two layers, mechanical stresses are induced in 
the tissue. Such a layered arrangement is commonly used to model various 
mechanical instabilities, e.g. the fragmentation of thin films |36|I37], and can be 
likened with the layered structure of plant leaves [TT] , with the rigid substrate 
representing the epidermis. 

The coupling between the two layers is described by a quadratic potential, 
which included into the energy (eq. 0|) yields 



E 



f i^f + f (e - e r {<i))f + M) + \{u~ 



dx. (9) 



The parameter k represents the strength of the coupling, and 77 is a measure of 
the mismatch between the two layers, defined by L s = (1 + r))L, where the L s 
and L are the lengths of the substrate and of the tissue layer in an unstressed 
state. 

To understand how the phase transition can be triggered, it is useful to 
compute the strain field when the tissue is in its initial state {<fi = 0) d. The 
equation of mechanical equilibrium (eq. [T]) yields 

Hd xx ii - k(u - r/x) = 0, (10) 

which is linear in u and admits solutions of the form 

u = r\x + Ae^ + Be~^ . (11) 

In this equation, A = \fjijk is a characteristic elastic length scale of the system, 
and A and B are constants that are determined by the boundary conditions. 
Here, the edges of the tissue are free of stresses, and d x u — at x = ±L S . 
Computing A and B and deriving eq. 1111 with respect to x, we find that the 
strain in the tissue is given by 

*=*( 1 - J *iPm)- (12) 

V cosh( J L s /A)/ 



3 Similar calculations could be carried out for the states of equilibrium in which the two 
phases coexist if the width of the interfaces is neglected. 
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Figure 3: A simulated pattern with three collapsed zones. The solid line and 
the dashed line represent the strain and the stress, respectively. The parameters 
of the model are fi = 1, e rl = -l,a= -0.25, (3 = 0.25, D = 4 x 10~ 4 , k = 16, 
L = 1, and r\ = —.5. With these values, the elastic length scale is A = .25, the 
instability threshold is e c o = —0.25 and the plateau stress is a p = 0. Accordingly 
the strains in the two phases when they coexist are equal to the rest strains in 
their rest configurations, i.e. e p o = and e p i = —1. 

If the system is large compared to the elastic length scale A, the strain is close 
to the strain imposed by the substrate everywhere except at the edges, and the 
stresses exhibit a broad plateau. If the system size is comparable to or smaller 
than the elastic length scale, the stresses have a more localized maximum at the 
center of the system, and the value of this maximum increases with the size of 
the system. 

One way instabilities can be induced is by increasing the mismatch between 
the two layers (this is analogous to the situation of a shrinking film on a sub- 
strate). From the above, it appears that instabilities can also be induced by 
growth with the mismatch between the two layers being fixed, that is as result 
of uniform growth. The latter regime will always be considered when simulating 
patterning in growing tissuefl 

We now turn to a numerical analysis of the evolution of the system. So 
far, we specified only its energy, which determines its states of equilibrium and 
their stability. To determine the new state of equilibrium reached after an 
instability occurs, its dynamics must be specified. A reasonable assumption is 
that mechanical relaxation is much faster than the transition between the two 
states of the tissue, so that at any time the system is at mechanical equilibrium: 

5E 

= — = -/J,d x {d x u - e r {4>)) + k(u - r/x) (13) 
ou 

Dynamics need thus be specified only for the phase field. We assume the latter 

4 If the two layers were instead assumed to grow at different rates, then the mismatch 
would become arbitrarily large with time. 
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obeys a relaxation equation of the form 

r) F 1 

d t <P~-—=DA(b~d (f> e-d (f> f. (14) 

At each time step of the overall evolution of the system, i.e. of its growth, this 
relaxation equation is solved assuming the system is at mechanical equilibrium, 
until an equilibrium is reached 0. And this state of equilibrium satisfies the 
equations^: 

f = fid x (d x u - e r {4>)) - k(u - -qx) . . 

\ = DA^-d^e-d^f. [L0) 

Most of the time, no instability occurs, and the new state of equilibrium is 
very close to the previous one. If an instability occurs, the phase field evolves to 
a new state that contains a different number of collapsed regions. Fig. [3] shows 
a typical equilibrium configuration. 

The results obtained in various conditions are shown in fig. 0] In the first 
example (fig. BJ a ))j the tissue has a fixed size and the mismatch between the 
two layers is progressively increased. A limited number of collapsed regions are 
formed and grow to occupy the entire tissue. In the second example (fig-Bib)), 
the system grows with the mismatch remaining constant (L = e~ Lq and L s = 
e~L s o). In that case, instabilities occur with time in an approximately self- 
similar cascade. Note that in both of these examples, the existing collapsed 
zones readjust after each new instability, and the system relaxes to a regular 
pattern that minimizes its energy. In the process, existing collapsed zones may 
change size or position, which is possible because the transition between the 
phases is reversible. In the case of a growing system (fig. 131b)), this means that 
the pattern does not retain any trace of the history of its formation. 

Now, if the transition to the collapsed state represents tissue differentiation, 
it should be irreversible, precluding the above rearrangements. A simple way 
of implementing irreversibility is to make the evolution of the phase field </> 
irreversible. Assuming that the transition must progress to a certain point 
before it becomes irreversible, we allow <f> to increase or decrease when it is 
smaller than a certain threshold (1/2), and only to increase when it is larger 
0. As shown in fig. [4jc), this yields hierarchical patterns, which reflect the 
history of the system. The collapsed zones grow along with the system, so that 
earlier-formed zones are larger than more recent ones. 

So far, we assumed that the tissue was perfectly homogeneous, and it may 
be of interest to examine how disorder, which is inherent in biological tissues, 

5 The equations of mechanical equilibrium (a set of linear equations) are solved at each 
step of the relaxation. The relaxation is stopped when the rate of change of the phase field 
falls below a certain threshold. 

6 We are thus assuming that there are three well separated time scales T mec h <C r^iff <C 
r growthi where T mec ^, Tdijj, and T growt h are the characteristic time scales corresponding to 
mechanical relaxation, cell differentiation, and growth, respectively. 

7 Such an ad hoc rule would clearly be unsatisfactory to describe a physical system, but 
in the case of biological systems, it can be seen as a necessary simplification of the complex 
phenomena involved. 
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would affect the outcome of the model. As described earlier, disorder can be 
incorporated as a multiplicative noise term in the energy of the system. In- 
stabilities then no longer occur exactly at the locations where the stresses are 
largest. Fig. QJd) shows an example of the resulting irregular patterns. Con- 
sistent with theoretical arguments and experimental results concerning crack 
patterns [37J IMj , we find that the effect of disorder depends crucially on how 
the distance between collapsed regions compares with the elastic length scale A. 
If this distance is much larger than A, then the stresses exhibit a broad plateau 
and the location of new collapsed zones is very variable, being governed essen- 
tially by the weak spots due to disorder. On the other hand, if the distance 
between collapsed regions is comparable to or smaller than A, then the stresses 
have a more localized maximum and the location of new collapsed zones is close 
to the location of this maximum. The latter regime is thus more favorable for 
obtaining regular and reproducible patterns. It is worth noting that in growing 
plant leaves, the distance between veins is of the same order as the thickness 
of the leaf [TOj (which sets the characteristic length scale for elastic coupling 
between veins). 



3 Two-dimensional model 

We now turn to the extension of the model to two-dimensional systems. In this 
context, the one-dimensional patterns described in the previous section can be 
interpreted as parallel stripes of localized deformation. The more general pat- 
terns that we expect to describe would take the form of networks of stripes of 
different orientations, corresponding to different principal directions of deforma- 
tion. It is thus clear that a scalar parameter is no longer sufficient to describe 
the collapsed state. Since the difference between the two phases lies only in their 
rest configurations, we choose here to identify the phase field tfi and e r , which in 
2D is a tensor describing the local rest configuration. Likewise, the mismatch 
between the two layers is now described by a tensor r\. The 2D analog of eq. [9] 
thus takes the form 



E 



y |Ve r | 2 + ^Tr(e - e r ) 2 + M |e - e r | 2 + f(e r ) + k(u - 



dS. (16) 



While r\ may be anisotropic, we assume that the tissue has no preferential di- 
rection. Accordingly, / can be expressed as a function of the invariants of e r , 
e.g. its trace and its norm. 

To obtain a behavior similar to that of the one-dimensional model, we con- 
sider potentials that have a minimum for e r = and a minimum of the same 
depth for values of e r corresponding to a smaller natural volume (i.e. Tre r < 0). 
Among all possible potentials satisfying this constraint, we have chosen the 
following three illustrative examples. In the first, the secondary minimum is 
reached for all deformations having a given amplitude (norm) , regardless of the 
shear. The second one specifically favors uniaxial deformations. In contrast, 
the third one has a secondary minimum for isotropic compression. 
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Figure 4: Patterns obtained in a one-dimensional system. These plots show 
the evolution of the phase field, with the collapsed regions appearing in black. 
Model parameters are given in fig. ([3]). (a) Increasing mismatch between the 
two layers, with the size of the tissue being fixed (L = 1). (b-d) Growing system 
with a fixed mismatch (77 = —0.275). The initial size of the tissue is Lg = 0.25. 
(b) Reversible evolution. Note that the collapsed regions appear to shrink only 
because the total size of the system increases, (c) Irreversible evolution, (d) 
Same as (c) with disorder. 
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Figure 5: Gray scale plots of the different potentials /(e r ) used. t\ and e 2 
are the eigenvalues of e r . Darker regions correspond to lower values, (a) The 
potential defined by eq. [XT] favors compressive strains of a given amplitude, 
regardless of the shear, (b) The potential defined by eq. [18] and the first set 
of coefficients given in the text favors uniaxial compression, (c) The potential 
defined by eq. [18] and the second set of coefficients favors isotropic compression. 



More specifically, the first potential is defined by 

/(e r ) = -(| ej .|~l/2) 2 + (| er |-l/2) 4 , (17) 



which has minima for e r = and for \e r \ = 1 (see fig. 5(a) I. The second and 



third examples are defined as polynomials of the invariants of e r , i.e. 

f(e r ) =^c ii (Tre r ) i |e r | 2 ^. (18) 

The two sets of coefficients used are c 2 o = .4, cqi = .2, C30 = .7, c\\ = .5, C40 = 
.4, C21 = .1, C02 = .1, which corresponds to minima for uniaxial deformations (see 
fig- |5(b)[ ) and c 20 = .25,c i = .5,c 30 = l,cn = 0, c 40 = -3 75 , c 2 i = 0, c 02 = .5, 
which leads to a minimum for isotropic compression (see fig. |5(c)[ ) Q 

In contrast with the one-dimensional system show in fig. [2] in which the 
edges were free of stresses, most of the 2D simulations were done with periodic 
boundary conditions, which allow the mechanical equilibrium equation to be 
efficiently solved in Fourier space. In this case, the phase field is initially chosen 
to be random, to break the symmetry of the system. The initial conditions 
(mismatch or size for a growing system) are chosen such than an instability 
immediately develops, preventing relaxation to a uniform state. 

As in the one-dimensional case, we first consider reversible evolution in a 
non-growing system (fig. [6] Here, the mismatch between the two layers is di- 
rectly set to its final value, and the system is allowed to relax to equilibrium. 
The first two potentials (shown in figs.[5](a) and (b)), which allow or favor uni- 
axial compression, yield very similar patterns of interconnected stripes, which 
could be likened to crack networks in directional growth systems such as basalt 



8 Note that the potential defined by eq. 1 171 cannot be written in this form. 
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Figure 6: Typical patterns observed in a layered two-dimensional system with 
the different potentials. Periodic boundary conditions were used and the initial 
state was randomly perturbed. Gray levels correspond to the trace of the rest 
configuration e r , with darker areas corresponding to collapsed regions. Poten- 
tials that allow or favor uniaxial compression (a and b, corresponding to the 
potentials of figs. 5(a) and 5(b) respectively) yield reticulate patterns, while 



the potential that favors isotropic compression produces discontinuous domains 
arranged in a hexagonal lattice (c, corresponding to the potential of fig. 5(c) I. 



columns |39j . In contrast, the potential shown in fig. 5(c) which favors isotropic 
compression, produces islands of localized deformation arranged in an hexago- 
nal lattice, which are similar to some structures observed in solid-solid phase 
transitions [57]. The characteristic length of the above patterns, as in ID, is set 
by the elastic length scale A = J fi/k. 

The difference between the two potentials that produce reticulate patterns 
becomes more apparent when considering a growing system (and a reversible 
evolution) , which makes it possible to observe individual nucleation events lead- 
ing to the formation of new collapsed regions. Indeed, in the case of the potential 



of fig. 5(a) new collapsed regions take the form of trijunctions that split an un- 
collapsed region into thre^l, while in the case of the potential of fig. 5(b) they 



occur as single stripes that divide an uncollapsed region across its longest ex- 
tension 0. This qualitative difference can be explained by the fact that the core 
of a trijunction is compressed in both directions, a state that is disfavored by 



the potential of fig. 5(b) 



Finally, we consider irreversible evolution in a growing system, which we 
found to lead to hierarchical patterns in ID. In that case, irreversibility was 
implemented by forcing the evolution of the scalar phase field to be monotonic 
once a certain threshold was exceeded. A straightforward extension to two di- 
mensions would be to apply the same rule to the trace of the rest configuration 
tensor e r . However, this would allow the principal directions of the rest config- 



9 A typical example in the case of an irreversible evolution can be seen in fig. ISja) 
10 A typical example can be seen in the case of an irreversible evolution can be seen in 

fig. He) 
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Figure 7: Patterns generated by irreversible evolution in growing systems with 
different potentials and boundary conditions (periodic boundary conditions or 
circular domains), (a) and (c) were obtained with the potential of fig. 5(a) (b) 
and (d) with the potential of fig. |5(b)| As expected, hierarchical structures are 
formed. However, they are disrupted by the tensions that accumulate within 
and around the collapsed regions as they grow (see text and fig. [5]) . Notice the 
lighter spots or stripes within the collapsed regions. 



uration in the collapsed phase to vary (and simulations indicate that this can 
indeed occur). Instead, we wish the establishment of these special directions to 
be definitive, as is the transition to the collapsed state (think of the elongation 
of vascular cells in plant leaves). Accordingly, we impose that, once the com- 
pression in some direction has exceeded a certain threshold, it can no longer 
decrease EE 

As in ID, we consider a growing system with a fixed mismatch, and find 
that hierarchical structures are obtained. As shown on fig. the shape of the 
patterns depends both on the potential used and on the shape of the growing 
domain (for the earliest structures that are formed). However, we find that the 
simple model presented here produces unexpected artifacts, such as regions of 
lower compression within the collapsed zones, and open-ended collapsed regions 
that fail to reconnect with previously formed structures. 

As illustrated by fig. these effects are caused by tensions that accumulate 
in the collapsed structures as they grow. When the transition to the collapsed 
state is reversible, the collapsed regions tend to an equilibrium width such that 
compressive stresses are relaxed in their neighborhood (i.e. for an ideal, in- 

11 To this end, we define a "maximum deformation tensor" e m that keeps track of the 
maximum past compression in any direction, i.e. the tensor e r (t) — e m (t') is positive for all 
times t < t'. From this maximum deformation tensor, we define an "irreversible deformation 
tensor" e;. In a basis that diagonalizes e m , ti = diag(/(Ai), /(A2)), where e m = diag(Ai, A2). 
The function / is defined by f(x) = 0.9x if x < —0.5 and f(x) = 00 if x > 0.5, so that 
compressions in excess of .5 are stored in e^. Irreversibility is implemented by imposing that 
the tensor e r — Ei is negative. The factor .9 in the definition of / was introduced to prevent 
numerical instability. 
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finitely long and straight stripe, the stress at the interface has a fixed value 
equivalent to the plateau stress a p of the one-dimensional model). In contrast, 
in the irreversible case, the width of each collapsed region grows along with the 
system and exceeds its equilibrium width. Instead of just releasing the com- 
pression in their vicinity, these growing collapsed regions give rise to tensile 
stresses. As can be seen on fig. [JJ these tensile stresses both destabilize the 
collapsed regions, and prevent the progression of new collapsed regions towards 
them. It must be noted that such tensile stresses are also generated in the ID 
model. However, they do not have such dramatic effects in that case, because 
new collapsed structures always form away from existing ones. 

4 Discussion 

In this article, we explored a possible mechanism for biological patterning by 
mechanical stresses. This mechanism, which involves the stress-mediated tran- 
sition of a tissue between two states, was shown to yield a response similar to 
that of a non-linear elastic material, justifying an analogy with mechanical in- 
stabilities. The formation of patterns was driven by coupling the tissue to a 
rigid substrate, which could represent a stiffer tissue such as the epidermis of 
plant leaves. As in the mechanical instabilities of layered physical systems, such 
as the wrinkling of a film bound to a substrate, this elastic coupling introduces 
a characteristic length scale that governs the size of the patterns. 

We first analyzed a one-dimensional system. In this case, the model involves 
a small number of parameters, and it is possible to describe the equilibrium pat- 
terns analytically. Numerical simulations showed that regular patterns, which 
are largely independent of the history of the system, are obtained when the 
transition between the two states of the tissue is reversible. In contrast, when 
this transition is made irreversible to represent tissue differentiation, the history 
of the system is retained in its final state. In a growing system, new structures 
keep forming over time, yielding hierarchical patterns. This is very similar to 
the development of leaf venation networks, in which veins of different orders 
form successively over the course of leaf growth |16| . 

In two dimensions, because of the tensorial nature of elastic fields, the poten- 
tial that characterizes the preferred states of the tissue can have many different 
forms. While we did not carry out a systematic analysis, the different exam- 
ples considered suggest that reticulate patterns can readily be obtained when 
the potential allows or favors uniaxial deformations of the tissue. In contrast, 
islands of deformed tissue are obtained with a potential that favors isotropic 
deformations. This suggests that a mechanical model of leaf venation pattern- 
ing could account for both wild-type patterns and the disconnected vasculature 
observed in certain mutants [40] , 

A very similar model of leaf vein patterning was recently proposed in [17] . 
In that approach, however, the different states of the tissue are characterized 
by a single scalar order parameter. The transition to the collapsed state is 
associated with an isotropic deformation, and the collapsed tissue is assumed 
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Figure 8: (Color online) Effect of irreversibility on the mechanical state of 
the system and pattern formation. All figures represent growing systems with 
periodic boundary conditions, (a) and (b) were obtained with the potential of 



fig. 5(a) (c) and (d) with the potential of fig. 5(b) Evolution was reversible in 



(a) and (c), irreversible in (b) and (d). Collapsed regions (defined by thresh- 
olding the value of Tr(e r )) are shown in blue. Gray levels indicate compressive 
stresses (lighter regions are more strongly compressed), and regions appearing 
in black are under tension. In the reversible case, the stresses decrease near the 
collapsed regions, yet generally remain compressive. In the irreversible case, this 
is also true for recently formed structures. In contrast, the larger, older struc- 
tures generate tensile stresses. These tensile stresses hinder the progression of 
newly formed structures, yielding disconnected patterns. 
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to have a reduced shear modulus to allow the formation of reticulate patterns. 
The elongation of vascular cells is described as a change of their rest volume 
combined with a strong elastic shear. One benefit of this approach is that it 
involves a much smaller number of parameters. However, it is unclear that 
it gives a realistic description of the mechanical state of vascular cells. While 
it comes at the cost of a greater complexity, we find it more appropriate to 
describe the elongation of vascular cells as a change in their rest shapes. This 
also makes it possible to address the effect of variations in these rest shapes 
on the overall pattern. For instance, as mentioned above, our results suggests 
that failure to elongate properly can lead to disconnected patterns such as are 
observed in some mutants. 

From the preceding results, we expected that the introduction of irreversibil- 
ity in a growing two-dimensional system would lead to hierarchical, reticulate 
patterns. However, we found that strong tensile stresses develop in collapsed 
regions as they grow beyond their initial, equilibrium size, and dramatically 
affect the evolution of the system. This is a consequence of the very simple 
growth laws assumed in our model. The rest configuration of the tissue changes 
when it switches between its two states, but growth in each of the two states 
is uniform and proceeds at a constant rate, allowing residual stresses to accu- 
mulate. In reality, the growth of biological tissues depends on the mechanical 
stresses to which they are subjected [3TJ 22] , which can be expected to limit the 
buildup of residual stresses. Incorporating this dependence would be essential 
to further investigate biological patterning by mechanical forces. Another per- 
spective would be to integrate the mechanical approach developed here with the 
biochemical factors of tissue differentiation. For instance, the differentiation of 
vascular cells in plant leaves involves the expression and polar localization of 
auxin carrier proteins [T3] , and it would be of interest to investigate how these 
processes are connected with the elongation of vascular cells. 
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